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Abstract 

In this paper we consider satellite trajectories in central force field with quadratic 
drag using two formalisms. The first using polar coordinates in which the angular 
momentum plays a dominant role. The second is in Levi-Civita coordinates in which 
the energy plays a central role. We then unify these two formalisms by introducing 
polar coordinates in Levi-Civita space and derive a new equation for satellite orbits in 
which energy and and angular momentum are on equal footing and thus characterize 
the orbit by its two invariants. In the second part of the paper we derive in Levi- 
Civita coordinates a linearized equation for the relative motion of two satellites whose 
trajectories are in the same plane. We carry out also a numerical verification of these 
equations. 


1 Introduction 


The accurate computation of satellites orbits has been the subject of numerous monographs 
[1-6] and research papers [7-9] (to name a few). In particular the effect of the Earth oblate¬ 
ness [6,17,18] and drag forces on satellite orbits [10,11,25] were the subject of some recent 
publications [19-21], 

Since satellite orbits in a central force field are in a plane it is natural to use polar 
representation for the orbit equations. In this formalism the angular momentum of the 
satellite emerges as a key variable for the derivation of the equations of motion. However 
in 1920 [6,23] Levi-Civita introduced another two dimensional formalism whose primary 
objective was to regularize the equations of motion near collision which is advantageous 
from computational point of view. To this end a coordinate transformation was defined and 
in the resulting equations of motion, energy emerges as the major quantity. A generalization 
of Levi-Civita coordinates to three dimensions was made by Kustaanheimo and Stcifel (KS- 
coordinates) [6,24,25] 

Related to the problem of satellite trajectory determination is the issue of relative motion 
of two satellites in orbit and the rendezvous problem. [12-16]. This problem was considered 
in several settings. In particular the rendezvous problem in the non-central force field of 
an oblate body was addressed in [18]. However as far as we know this problem was not 
addressed even in two dimensions using Levi-Civita coordinates which mitigate to some 
extent the singular nature of these equations. 

Our primary objective in this paper is to introduce polar coordinates in Levi-Civita plane 
and derive a new equation which charaterize the motion of a particle in a central force field 
in terms of its two invariants viz. energy and angular momentum. These invariants appear 
on equal footings in this equation. Thus this equation might find applications in classical 
mechanics and in satellite theory, mission planing and control. Our secondry objective in this 
paper, is to recast in Levi-Civita coordinates the relative motion of two satellites orbiting in 
the same plane thereby mitigating the sigular nature of these equation when these satelittes 
are in close proximity. 

Within this framework we consider also the dissipative effects of quadratic drag on the 
the angular momentum and energy of a satellite and their impact on it orbits. This problem 
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was addressed in various context by several authors [16,18,19,26,27] in the past. 

The plan of the paper is as follows: In Sec 2 we review the general theory for satellite 
trajectories in a central force held and derive the orbit equation under the action of quadratic 
drag. Sec 3 provides a short review of Levi-Civita formalism. In Sec 4 we introduce polar 
coordinates in Levi-Civita plane and derive the orbit equation in these coordinates. In Sec 
5 we derive a linearized equation for the relative motion of two satellites in the same plane. 
Sec 6 carries out some numerical simulations to verify the accuracy of the formulas that were 
derived in Secs 4, 5. We end up with some conclusions in Sec 7. 

2 Angular Momentum Representation of the Orbit 

In this paper we consider satellites in a central force held with quadratic drag. The general 
equation for the orbit of a satellite under these assumptions is 

R = —f(R)K - g(a, R)(R • R) 1/2 R, R — |R| (2.1) 

In this equation R is the radius vector of the satellite from the center of attraction and a is 
a parameter that lumps together drag constant and the atmospheric density proportionality 
constants. Differentiation with respect to time is denoted by a dot. We assume that / and 
g are differentiable on the domain of R under consideration. 

Taking the vector product of (12.ip with R on the left and introducing the angular mo¬ 
mentum vector 

J = R x R (2.2) 

we obtain 

j = -g(a,R)(R- R) 1/2 J (2.3) 

It follows from this equation that J is always parallel to J and therefore J has a hxed 
direction. As a consequence the motion is in a hxed plane which we take, without loss of 
generality, to be the x — y plane. Introducing polar coordinates R and 6 in this plane the 
angular momentum vector can be written as 

J = R 2 9ej (2.4) 
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where ej is a unit vector in the direction of J. Hence (12.3p can be rewritten as 

j 


J 


= -g(a,R)( RR) 1/2 


(2.5) 


where J = |J|. Substituting this result in (12. ip and dividing by J yields 

d /r\ f(R )R 


dt \ J 


+ 


J 


= 0 


i.e 


R= _j rim* 


J 


In polar coordinates (12. ip becomes 


R6 + 2R9 = —g(a, R)(R 2 + R 2 6 2 ) 1/2 R9 

R - R9 2 = —f(R)R - g(a, R)(R 2 + R 2 9 2 ) 1/2 R 
Dividing the first equation by R9 and integrating we obtain 

J = R 2 9 = hex p j g(a, R)(R 2 + R 2 9 2 ) 1 ^ 2 dt 


( 2 . 6 ) 

(2.7) 

( 2 . 8 ) 
(2.9) 

( 2 . 10 ) 


where h is an integration constant. 

Using (12.4p to change the independent variable from t to 9 we obtain after some algebra 
the orbit equation 

R" f R '\ 2 f(R)R 4 

lf- 2 U) + %V = 1 (2 ' n) 

where primes denote differentiation with respect to 9 

Equation (12.lip is the orbit equation for the motion of the satellite in the ”angular 
momentum representation”. 

When f(R) R = VD (12.ip becomes 


R = -W - g(a, R)( R • R) 1/2 R 


Taking the scalar product of this equation by R we obtain 


( 2 . 12 ) 


(R.R) = —(W,R) -g(a,R)(R- R) 3/2 . 


(2.13) 
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This can be rewritten as 


(2.14) 


d£ 

dt 


—g(a, i?)(R • R) 3/2 . 


where E is the particle energy 


E 


V + 


(R,R) 

2 


Eq. (I2.14p gives the rate of the particle energy decay due to the dissipative effects of the 
drag force. 


3 Energy Representation of the Orbit 

Energy is an invariant which characterize satellite motion when no drag is present. To 
take advantage of this fact Levi-Civita[ ] introduced a two dimensional formalism (which 
was latter generalized to three dimensions[]) in which Energy plays a central role (and has 
additional advantages when collisions are present). 

We present here a short summary of this formalism [ j 


3.1 Levi-Civita formalism 


There are in the literature excellent expositions of Levi-Civita formalism [6]. Here we present 
a short overview of this formalism for completeness. 

To begin with introduce a ’’fictitious time” s which is defined by the relation 


ds dt 


(3.1) 


It is then easy to see that for R = (x, y ) 


R = - ^R'R'- (3-2) 

where primes denote differentiation with respect to s. Furthermore the velocity v satisfies, 

u 2 = (R,R) = -^(R',R') (3.3) 

Using (13.21) and (13.31) the equation of motion (12.11) becomes 

“ lv R ' R ' = ~ ~ 7 ^ R/ ) V2 R/ ( 3 - 4 ) 
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For a particle of unit mass whose equation of motion is (12. ip where f(R )R = W and 
g(a, R) = 0 (no drag) the energy E is conserved and 

E = V+ V ^ = (3.5) 

When this particle is in the gravitational field of a point mass M, 


m 



h 

R 


(3.6) 


and 

E — t k + 27p ( r ' r '>- < 3 ' 7 > 

where fi = GM and G is the gravitational constant. The rate of change in E in this 
coordinate system can be obtained from (12. 14[) by a change of variables 


dE 

ds 


R 2 


(R/, R') 3/2 


(3.8) 


Next define a transformation from the (x,y) coordinate system to a new one (tq,-u 2 ) (Levi- 
Civita coordinates) which is defined by the following relations, 


x = u\ — u\i y = 2 u\U 2 - 


(3.9) 


(We shall refer to this as the U-plane). 

Introducing the scalar product of any two vectors wy, w 2 as 

(wi, w 2 ) = w / 1 w 2 . 


it follows that R = (u, u) = |u| 2 where u = (iq, u 2 ). Next we introduce Levi-Civita matrix, 


L(u) = 



-u 2 

Ml 


Observe that the transpose L T and the inverse L 1 
relationships 

L t (u)L(u) = RL, L~\ u) = 


of this matrix satisfy the following 
■ ii T (u) (3.10) 
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where I is the unit two dimensional matrix. Moreover for any two vectors u, v we have the 
following identities 


L(u)v = L(v)u, (u, u)L(v)v + (v, v)L(u)u = 2(u, v)L(u)v (3.11) 

and 

L(u)' = L(u') (3.12) 

It is then easy to see that for R 

R = L(u)u, R' = 2L(u)u', R" = 2L(u)u" + 2L(u)'u' = 2L(u)u" + 2L(u / )u / (3.13) 


To convert (j3.4[) to an equation in U space we use (13.131) and (13.111) with the vectors u, u'. 
After some algebra we obtain 

u" - ^ - U ) u = {-f(R)n - kg{oi, R)( u, u)” 3/2 (u', u') 1 / 2 u'} (3.14) 

(u, u) 2 k J 

where R has to be replaced by (u, u). When f(R) is given by (13.61) this equation becomes 


u" + ’ u ) u _ _2 g( a , R)( u, u) 1 ^ 2 (u / , u') 1 / 2 ^ 


When no drag forces are present the particle energy eq. (13.71) becomes 

/.t — 2(u', u') 


E = 


(u,u) 


Using this expression for E in (I3.15P we have 


u" — \^Eu = —2 g(a, R)( u, u) 1 / 2 (u', u') 1 / 2 ^. 


(3.15) 


(3.16) 


(3.17) 


However observe that due to the dissipative nature of the drag force, E is not constant under 
the present settings. 


4 Polar Representation of the Orbit Equation 

Polar coordinates representation of satellite orbit equations in the x — y plane offers several 
advantages over their Cartesian counterparts. Motivated by this observation we develop in 
this section a representation of (I3.15P using polar coordinates using Levi-Civita coordinates. 
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4.1 Polar geometry in the U-plane 

We introduce polar coordinates (u, 0) in the U-plane in a manner similar to polar coordinates 
in the the x — y plane. 

u = (u, u) 1//2 , (j) — tan~ 1 — . (4.1) 

U\ 

The relationship between these variables and (R , 9) is given by 

R = (u, u) = w 2 , 9 = 2<j). (4.2) 

Furthermore in parallel to the definitions of the radial and tangential unit vectors in x — y 
plane 

e r = (cos 9, sin 0), eg = (— sin#, cos 9), (4.3) 

we define in the U-plane 


e u = (cos 0, sin 0), = (—sin 0, cos 0). 


(4.4) 


Using M we find that 


e r = cos 0 e u + sin 0 e^, eg = — sin 0 e u + cos 0 e«y 
Since u = ue u we have the following formulas for the derivatives of u, 

u' = u'e u + ri0'e^,, u" = {u" - u((j)') 2 )e u + {u<j>" + 2u'4>')e ( j > . 


(4.5) 


(4.6) 


4.2 Derivation of the New Orbit Equation 

Using (14.6p the orbit equation (j3. 14(1 becomes 

| u" - n(0 / ) 2 ]e. u + [u(j>" + 2u'0 / ]e ?i - 


(u 7 , u') 


— 


U 


(u,u ) 2 


{-f(R)ue u - 4g(a,R)(u,u) 3/2 (u',u ') l/2 [u'e u + u0'e^]} 


This yields the following two equations for the tangential and radial components 

ucj)" + 2 u'cj)' = —2 g(a, R)u 2 { u', u , ) 1//2 0 / , 


(4.7) 


(4.8) 


u" - n(0') 2 - 


fu', u') 


f{R)u 5 


u 


2 


2g(a, R)u( u\ u' l IJ (/. 


(4.9) 












where 


(u', u') = ( u') 2 + u 2 (0') 2 . 


Multiplying (14.81) by u, dividing by u 2 (j)' and integrating we obtain 

L — u 2 <p' — hi exp (—2 / g(a, u)u(u\ u') 1 / 2 ds 


(4.10) 


(4.11) 


Hence 

d L d 

ds v? d(j) 

We note that L defined in (14.10j) is equal to the angular momentum J up to a constant. 
In fact using (14. 2 p we have 


2jJ _ VO 


L = u z (p' = u z R -- = 

r dt 2 

Using (14.lip to change the variable from s to 0 in (14.91) we obtain after a long algebra the 
following orbit equation, 


1 d 2 u 2 ddu\ 2 -u 2 (u', u') 1 / 2 f{u)u 8 

u~d4?~ v?\d4>) L 2 ““ 2 L 2 

When f(R) is given by (13.61) this becomes 

1 d 2 u 2 /du\ 2 u 2 [/i — 2(u', u') 1 / 2 ] 

udcf) 2 u 2 \d(p) 2 L 2 

Using (13.15ft this equation becomes 

1 d 2 u 2 Z'du\ 2 Eu 4 

u # 2 _ n 2 V#/ _ 2L^ “ 


(4.12) 


(4.13) 


(4.14) 


When there are no drag forces E, L are constants and eq. (14.141) may be referred to as 
the ’’Energy Angular Momentum” orbit equation in Levi-Civita coordinates. In absence of 
dissipative forces this equation characterize the orbit by its generic invariants. 


5 Relative motion of Satellites 

In this section we derive the equations for the relative motion of two satellites in orbit in a 
central force field. 
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5.1 Linearized Equations in Physical Space 

If the positions of the two satellites are denoted by Ri,R 2 then their respective equations 
of motion in a (conservative) central force field are 

Ri = -W(Ri), R 2 = -W(R 2 ) (5.1) 


where the dots represent differentiation with respect to time. The relative position of the 
second satellite with respect to the first is w = R 2 — Rx . This leads to the equation 

i?x + w = — VU(R 2 ). (5.2) 

Using (15.ip this can be rewritten as 

w = W(Ri) - W(R 2 ). (5.3) 

Assuming that |w| «Ki we can approximate 

VU(Ri) - VU(R 2 ) = V [V (Ri) - U(Ri + w)] 


by a first order Taylor polynomial in w. This leads to the following linear relative motion 
of the second satellite with respect to first in the inertial coordinate system attached to the 
central body center, 

w = — V(W(i?i) • w) (5.4) 


In particular if the motion is around a spherical body where V is given by (I3.6j) we have 

u 3u(Ri • wl 

W = -Jr +J V ^ 1 = F ' (5 ' 5) 

In a coordinate system rotating with the first satellite the relative-motion equation (15.411 
becomes [21] 

w + 2f2 xw + f2x(f2xw) + f2xw = F. (5-6) 


where f2 is the orbital angular velocity of the first satellite. The reduction of this formula to 
a system of ordinary differential equations for the motion of two satellites around an oblate 
body was carried in [18]. 
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We now consider this equation in the special case where the two satellites are in the same 
x — y plane. In this case 

f2 = (0,0,0), w = (w 1 ,w 2 ,0). (5.7) 

We have 

12 x w = Q(—w 2 , wi, Of, ft x (fi x w) = — 0 2 {wi, w 2 , 0 f, ft x w = 0(—w 2 , w\, 0) T 
(In the following we suppress the third component of the vectors). 


5.2 Relative Equation of Motion in Levi-Civita Coordinates 

We now introduce Levi-Civita transformation 


Wi — vf- v\, w 2 = 2viv 2 , r 2 = w\ + wl = (v, v) 2 , ~r =r^ 


(5.8) 


ds dt 

Due to the appearance of the vector (—w 2 , w \) in the equation of motion (15. 6 p we intro¬ 
duce 


L(u) = 


— U 2 U\ 


Ui —u 2 


We then have 


L(u)u = (-u 2 ,uif. 

The matrix Z(u) has the following properties 

L(u) T Z(u) = (u, u)J, L(u) t L(u) = (u, u)T, L^ 1 (u) t L(u) 


(5.9) 


u, u 


-LiufLiu) = T, 


(5.10) 


where 


r = 


o -1 

1 o 


Observe also that 


-Wo 


w'i 


= 2 L(v)V. 
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After some algebra similar to the the one in Sec 4 we obtain the following representation of 


eq. (|576|) 


Using 


(vVQ | r 2 9 2 


v, v 


■ . r 2 ti r 2 L T (v)F 

v + 2r6Tv H-Tv =-. 

2 2 


(5.11) 


« = i 


8" - ~(v,V)8' 

r 


eq (15. 11 j) becomes 


v"- 


v',v') , (O') 2 1 


(v,v) 


+ 


9 " - 


V, V 


r \ v + 2(6»)Tv = r 2 


L T (v)F 


(5.12) 


In particular when F is given by (15. 5 p the right hand side of this equation becomes 


2 Z/(v) F ar 6 ?>ur- T 

r—= -#=rV + • w)L t (v)R! 


2R\ 2R\ 


6 Numerical Verification 


6.1 Motion of a Satellite in an Exponential Atmosphere 


A common model for the earth atmosphere density p with height is 

p = Uiexp (^ R ° h R ^j (6.1) 

where C±, Ro, H are constants. For this atmospheric model 

. , f R() — R\ ( lin — U 2 \ . , 

g(a, R) = a exp f ———J = a exp i———\ (6.2) 

where the constant C\ was lumped with the drag coefficient a. Eq. (13.81) for the energy 
becomes 


ds 


^( u '. u ')3/2 

U 


exp 



(6.3) 


Similarly for L we have 


1 dL 
L ds 


—2au(u' ■ u') 1 / 2 exp 



(6.4) 
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Using (14. lip to change variables from s to <p in (I6.3[) . (16. 4 1) yields, 


dE 

dcp 

dL 

dcf) 


8 aL 2 (du du\ 


u° 


3/2 


Wm>) exp 


Uq — u 

H 


= —2 auL 


du du\ 
d(j) ’ dtp) 


1/2 


exp 


Un 


U 


H 


(6.5) 


( 6 . 6 ) 


The system of equations (14.14p . (16.5(1 . (16.6|) comprise of four equation in four unknowns and 
can be solved by numerical methods. Fig 1 is a plot of the solution of this system with 
Rq = 7000 km, 9 q = . HR (initial circular orbit), a = 3.10 10 , H = 88.667 and step error of 

V -^o 

10 -12 . On the same figure we plotted also the numerical solution of the system (12.8ft - (12.9p 
with the same parameters. We note that these curve are almost inditinguishable. Fig 2 
displays the difference between these curves which over ten periods remains less than 0.75m 
This difference is most probably due to the cumulative error in the numerical integration. 

To verify numerically the formula for the relative motion of satellites which was derived 
in the previous section we considered two satellites in circular orbit whose position at time 
t = 0 (in polar coordinates) is (i?i,0) and (i?2,0) with R± = 7000 km and R 2 = 6999 km. 
The angular velocities of these satellites respectively are 


OJi 

Hence their distance d at time t satisfies 



!±_ 

nr 


i = 1,2 


d 2 = R\ + Rl — 2RiR 2 cos(u;i — u 2 )t 

Fig 3 is a plot of the difference between this analytical expression for the distance and the 
numerical value obtained from the linearized formula (15.111) as a function of time. 


7 Conclusion 

In the first part of the paper we developed a new representation for the orbit equation of a 
satellite in terms of it natural invariants. In the second part a formula was derived for the 
relative motion of two satellites moving in the same plane. This formula can be generalized 
to the case where drag effects have to be taken into account and the orbits of the satellites 
are not in the same plane (using KS-formalism). 
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R as a function of 6 



Figure 1: : Illustrative trajectory for a satellite orbit using eq. (14. 14ft (red line) which is 
indistinguishable from the one obtained from fl2.8jl - fl2.9jl . 
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Difference in R as a function of 0 



Figure 2: : Diffrence between the trajectories in Fig 1 
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t (sec) 


Figure 3: : Diffrence between the analytic and numerical value of the distance between two 
satellites 
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